####################################
# Placebo RDROBUST
####################################

rm(list=ls())

#sink("~/Dropbox/Crime Chile/11_replication/08_rdrobust_placebo_crime.txt")

library(Hmisc)
library(ggplot2)
library(stargazer)
library(foreign)
library(rdrobust)
library(rdd)
library(readstata13)

################
# Prepare data 
################

# read data
d = read.dta13("~/Dropbox/Crime Chile/11_replication/local_crime_data_chile_2022january.dta")   
names(d)

############################
# RDrobust MSE
############################

# disturbing the peace (DTP)
summary(rdrobust(d$truidos_molestos_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA3a.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$truidos_molestos_s, x = d$margin,  h= 0.148, nbins = 100, subset = -0.148 <= d$margin & d$margin <= 0.148, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-1, 1), title = "                                       Robust CI: [-0.332 , 0.076]",
        y.label = "Disturbing the peace", x.label = "Margin of victory"))
dev.off()

# domestic violence against the elderly (DVE)
summary(rdrobust(d$tviolencia_adultomayor_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA3b.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$tviolencia_adultomayor_s, x = d$margin,  h= 0.125, nbins = 100, subset = -0.125 <= d$margin & d$margin <= 0.125, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-1, 1), title = "                                       Robust CI: [-0.083 , 0.389]",
        y.label = "Domestic violence against the elderly", x.label = "Margin of victory"))
dev.off()

# domestic violence against men (DVM)
summary(rdrobust(d$tviolencia_hombre_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA3c.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$tviolencia_hombre_s, x = d$margin,  h= 0.152, nbins = 100, subset = -0.152 <= d$margin & d$margin <= 0.152, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-1, 1), title = "                                       Robust CI: [-0.127 , 0.330]",
        y.label = "Domestic violence against men", x.label = "Margin of victory"))
dev.off()

# domestic violence against women (DVW)
summary(rdrobust(d$tviolencia_mujer_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA3d.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$tviolencia_mujer_s, x = d$margin,  h= 0.238, nbins = 100, subset = -0.238 <= d$margin & d$margin <= 0.238, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-1, 1), title = "                                       Robust CI: [-0.262 , 0.240]]",
        y.label = "Domestic violence against women", x.label = "Margin of victory"))
dev.off()

# domestic violence against children (DVC)
summary(rdrobust(d$tviolencia_ninos_s,d$margin,cluster = d$cluster,all=TRUE))

cairo_pdf(file="~/Dropbox/Crime Chile/11_replication/figureA3e.pdf", 
          width=7, 
          height=7)
(rdplot(y = d$tviolencia_ninos_s, x = d$margin,  h= 0.219, nbins = 100, subset = -0.219 <= d$margin & d$margin <= 0.219, 
        binselect="esmv", kernel="triangular", p=1, col.lines = "black", col.dots = "black", y.lim = c(-1, 1), title = "                                       Robust CI: [-0.334 , 0.155]",
        y.label = "Domestic violence against children", x.label = "Margin of victory"))
dev.off()

#############################
# Table A15 (first five rows)
#############################

summary(rdrobust(d$truidos_molestos_s,d$margin,cluster = d$cluster,all=TRUE))
pe1 = as.numeric(rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv1 = as.numeric(rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci1 = as.numeric(rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss1 = as.numeric(rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess1 = as.numeric(rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw1 = as.numeric(rdrobust(d$truidos_molestos_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

summary(rdrobust(d$tviolencia_adultomayor_s,d$margin,cluster = d$cluster,all=TRUE))
pe2 = as.numeric(rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv2 = as.numeric(rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci2 = as.numeric(rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss2 = as.numeric(rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess2 = as.numeric(rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw2 = as.numeric(rdrobust(d$tviolencia_adultomayor_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

summary(rdrobust(d$tviolencia_hombre_s,d$margin,cluster = d$cluster,all=TRUE))
pe3= as.numeric(rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv3 = as.numeric(rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci3 = as.numeric(rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss3 = as.numeric(rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess3 = as.numeric(rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw3 = as.numeric(rdrobust(d$tviolencia_hombre_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

summary(rdrobust(d$tviolencia_mujer_s,d$margin,cluster = d$cluster,all=TRUE))
pe4 = as.numeric(rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv4 = as.numeric(rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci4 = as.numeric(rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss4 = as.numeric(rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess4 = as.numeric(rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw4 = as.numeric(rdrobust(d$tviolencia_mujer_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

summary(rdrobust(d$tviolencia_ninos_s,d$margin,cluster = d$cluster,all=TRUE))
pe5 = as.numeric(rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$Estimate[1])
pv5 = as.numeric(rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$pv[3])
ci5 = as.numeric(rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$ci[3,])
oss5 = as.numeric(rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$N[1] + rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$N[2])
ess5 = as.numeric(rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[1] + rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$N_h[2])
bw5 = as.numeric(rdrobust(d$tviolencia_ninos_s,d$margin,all=TRUE,cluster = d$cluster)$bws[1,1])

pe = c(pe1,pe2,pe3,pe4,pe5)
pv = c(pv1,pv2,pv3,pv4,pv5)
ci_lower = c(ci1[1],ci2[1],ci3[1],ci4[1],ci5[1])
ci_upper = c(ci1[2],ci2[2],ci3[2],ci4[2],ci5[2])
oss = c(oss1,oss2,oss3,oss4,oss5)
ess = c(ess1,ess2,ess3,ess4,ess5)
bw = c(bw1,bw2,bw3,bw4,bw5)
variable = c("DTP","DVE","DVM","DVW","DVC")
data = data.frame(variable,pe,pv,ci_lower,ci_upper,oss,ess,bw)

colnames(data) = c("","Point Estimate","Robust P-value","Robust 95% Confidence Interval lower bound","Robust 95% Confidence Interval upper bound","Overall sample size","Effective sample size","MSE bandwidth")
data

stargazer(data, summary=FALSE, rownames=FALSE, out="~/Dropbox/Crime Chile/11_replication/tableA15_part1.html")

sink()

